home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / SOR.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  1KB  |  49 lines
  1. PROCEDURE sor(a,b,c,d,e,f: gljmax; VAR u: gljmax;
  2.          jmax: integer; rjac: double);
  3. (* Programs using routine SOR must define the type
  4. TYPE
  5.    gljmax = ARRAY [1..jmax,1..jmax] OF double;
  6. in the main routine. *)
  7. LABEL 99;
  8. CONST
  9.    maxits=1000;
  10.    eps=1.0e-5;
  11.    zero=0.0;
  12.    half=0.5;
  13.    qtr=0.25;
  14.    one=1.0;
  15. VAR
  16.    n,l,j: integer;
  17.    resid,omega,anormf,anorm: double;
  18. BEGIN
  19.    anormf := zero;
  20.    FOR j := 2 TO jmax-1 DO BEGIN
  21.       FOR l := 2 TO jmax-1 DO BEGIN
  22.          anormf := anormf+abs(f[j,l])
  23.       END
  24.    END;
  25.    omega := one;
  26.    FOR n := 1 TO maxits DO BEGIN
  27.       anorm := zero;
  28.       FOR j := 2 TO (jmax-1) DO BEGIN
  29.          FOR l := 2 TO (jmax-1) DO BEGIN
  30.             IF (((j+l) MOD 2) = (n MOD 2)) THEN BEGIN
  31.                resid := a[j,l]*u[j+1,l]+b[j,l]*u[j-1,l]
  32.                   +c[j,l]*u[j,l+1]+d[j,l]*u[j,l-1]
  33.                   +e[j,l]*u[j,l]-f[j,l];
  34.                anorm := anorm+abs(resid);
  35.                u[j,l] := u[j,l]-omega*resid/e[j,l]
  36.             END
  37.          END
  38.       END;
  39.       IF (n = 1) THEN BEGIN
  40.          omega := one/(one-half*sqr(rjac))
  41.       END ELSE BEGIN
  42.          omega := one/(one-qtr*sqr(rjac)*omega)
  43.       END;
  44.       IF ((n > 1) AND (anorm < (eps*anormf))) THEN GOTO 99
  45.    END;
  46.    writeln('pause in routine SOR');
  47.    writeln('too many iterations'); readln;
  48. 99:   END;
  49.